library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
## 
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
#library(corrplot)
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)

RRPLOTS and flchain

odata <- flchain
odata$chapter <- NULL
table(odata$death)
## 
##    0    1 
## 5705 2169
rownames(odata) <- c(1:nrow(odata))
data <- as.data.frame(model.matrix(Surv(futime,death)~.*.,odata))

data$`(Intercept)` <- NULL
table(odata[rownames(data),"death"])
## 
##    0    1 
## 4562 1962
dataFL <- cbind(time=odata[rownames(data),"futime"],status=odata[rownames(data),"death"],data)
dataFL$time <- dataFL$time/365
colnames(dataFL) <-str_replace_all(colnames(dataFL)," ","")
colnames(dataFL) <-str_replace_all(colnames(dataFL),"\\.","_")
colnames(dataFL) <-str_replace_all(colnames(dataFL),":","_")
colnames(dataFL) <-str_replace_all(colnames(dataFL),"-","_")
colnames(dataFL) <-str_replace_all(colnames(dataFL),">","_")

trainsamples <- sample(nrow(dataFL),2000)
dataFLTrain <- dataFL[trainsamples,]
dataFLTest <- dataFL[-trainsamples,]

pander::pander(table(dataFLTrain$status))
0 1
1395 605 
pander::pander(table(dataFLTest$status))
0 1
3167 1357

Modeling

ml <- BSWiMS.model(Surv(time,status)~1,data=dataFLTrain,loops=1)
sm <- summary(ml)
pander::pander(sm$coefficients)
Table continues below
  Estimate lower HR upper u.Accuracy
age 0.05515 1.051 1.057 1.062 0.707
age_creatinine 0.05102 1.047 1.052 1.057 0.697
sample_yr_creatinine -0.00187 0.998 0.998 0.998 0.556
age_kappa -0.00748 0.992 0.993 0.993 0.712
age_flc_grp -0.00145 0.998 0.999 0.999 0.660
kappa 0.86847 1.937 2.383 2.932 0.656
flc_grp 0.09260 1.071 1.097 1.123 0.619
age_sexM 0.00199 1.001 1.002 1.003 0.528
kappa_lambda -0.02023 0.971 0.980 0.989 0.711
sexM_lambda 0.24150 1.133 1.273 1.431 0.544
lambda_flc_grp 0.03684 1.015 1.038 1.061 0.662
age_lambda -0.00363 0.994 0.996 0.999 0.722
sexM_kappa -0.19463 0.696 0.823 0.974 0.550
Table continues below
  r.Accuracy full.Accuracy u.AUC r.AUC
age 0.758 0.716 0.720 0.739
age_creatinine 0.742 0.711 0.682 0.737
sample_yr_creatinine 0.742 0.711 0.559 0.742
age_kappa 0.722 0.711 0.690 0.739
age_flc_grp 0.719 0.716 0.672 0.743
kappa 0.724 0.711 0.637 0.740
flc_grp 0.719 0.716 0.642 0.744
age_sexM 0.713 0.716 0.508 0.740
kappa_lambda 0.701 0.711 0.635 0.727
sexM_lambda 0.701 0.711 0.519 0.724
lambda_flc_grp 0.709 0.716 0.638 0.740
age_lambda 0.711 0.716 0.698 0.742
sexM_kappa 0.707 0.711 0.521 0.731
Table continues below
  full.AUC IDI NRI z.IDI z.NRI
age 0.744 0.08363 0.6485 15.80 15.428
age_creatinine 0.736 0.07762 0.7074 14.92 17.076
sample_yr_creatinine 0.736 0.06594 0.6725 13.09 15.899
age_kappa 0.736 0.00943 0.5048 8.70 11.434
age_flc_grp 0.744 0.00744 0.3111 8.45 6.887
kappa 0.736 0.01415 0.5305 8.24 11.904
flc_grp 0.744 0.00570 0.3186 7.64 7.088
age_sexM 0.744 0.00644 0.0392 4.59 0.808
kappa_lambda 0.736 0.00616 0.1729 4.18 3.565
sexM_lambda 0.736 0.00604 0.2266 4.05 4.770
lambda_flc_grp 0.744 0.00439 0.2409 3.22 5.039
age_lambda 0.744 0.00309 0.2833 3.15 5.911
sexM_kappa 0.736 0.00221 0.1881 2.27 3.881
  Delta.AUC Frequency
age 0.005213 1
age_creatinine -0.001161 1
sample_yr_creatinine -0.005483 1
age_kappa -0.002847 1
age_flc_grp 0.001594 1
kappa -0.003344 1
flc_grp 0.000299 1
age_sexM 0.003664 1
kappa_lambda 0.009041 1
sexM_lambda 0.011849 1
lambda_flc_grp 0.004082 1
age_lambda 0.002180 1
sexM_kappa 0.004739 1

Cox Model Performance

Here we evaluate the model using the RRPlot() function.

The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

timeinterval <- mean(subset(dataFLTrain,status==1)$time)

h0 <- sum(dataFLTrain$status & dataFLTrain$time <= timeinterval)
h0 <- h0/sum((dataFLTrain$time > timeinterval) | (dataFLTrain$status==1))

pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")
Initial Parameters
h0 timeinterval
0.161 5.73 
index <- predict(ml,dataFLTrain)
rdata <- cbind(dataFLTrain$status,ppoisGzero(index,h0))

rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataFLTrain$time,
                     title="Raw Train: FLC",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

As we can see the Observed probability as well as the Time vs. Events are not calibrated.

Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$OERatio),caption="O/E Ratio")

Quitting from lines 105-118 [unnamed-chunk-5] (flchain.Rmd) Error in [.data.frame: ! undefined columns selected Backtrace: 1. pander::pander(t(rrAnalysisTrain\(OERatio), caption = "O/E Ratio") 2. pander:::pander.htest(t(rrAnalysisTrain\)OERatio), caption = “O/E Ratio”) 3. pander::pandoc.table(res, caption = caption, …) 5. pander::pandoc.table.return(…) 7. base::sapply(1:ncol(t), function(x) is.numeric(t[, x])) 8. base::lapply(X = X, FUN = FUN, …) 9. pander (local) FUN(X[[i]], …) 11. base::[.data.frame(t, , x) Warning messages: 1: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 2: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 3: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 4: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 5: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 6: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 7: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 8: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 9: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter

pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
0.679 0.678 0.674 0.683 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
1.13 1.13 1.13 1.13 
pander::pander(t(rrAnalysisTrain$c.index$cstatCI),caption="C. Index")
C. Index
mean.C Index median lower upper
0.777 0.777 0.758 0.795 
#pander::pander(rrAnalysisTrain$c.index,caption="C. Index")
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.824 0.804 0.845 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.562 0.521 0.602 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.882 0.914 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.302 0.215 
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")

Quitting from lines 105-118 [unnamed-chunk-5] (flchain.Rmd) Error in t.default(): ! argument is not a matrix Backtrace: 1. pander::pander(t(rrAnalysisTrain\(RR_atP), caption = "Risk Ratio") 3. base::t.default(rrAnalysisTrain\)RR_atP)

pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 694.989012 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1305 189 436.0 139.89 505.22
class=1 214 76 65.9 1.54 1.73
class=2 481 340 103.1 544.04 664.87

Cox Calibration

op <- par(no.readonly = TRUE)


calprob <- CoxRiskCalibration(ml,dataFLTrain,"status","time")

( 13.74924 , 9120.585 , 1.200408 , 604 , 648.2787 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")
h0 Gain DeltaTime
0.267 0.758 13.7

The RRplot() of the calibrated model

h0 <- calprob$h0
timeinterval <- calprob$timeInterval;

rdata <- cbind(dataFLTrain$status,calprob$prob)


rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataFLTrain$time,
                     title="Calibrated Train: FLC",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

Calibrated Train Performance

pander::pander(t(rrAnalysisTrain$OERatio),caption="O/E Ratio")

Quitting from lines 156-169 [unnamed-chunk-8] (flchain.Rmd) Error in [.data.frame: ! undefined columns selected Backtrace: 1. pander::pander(t(rrAnalysisTrain\(OERatio), caption = "O/E Ratio") 2. pander:::pander.htest(t(rrAnalysisTrain\)OERatio), caption = “O/E Ratio”) 3. pander::pandoc.table(res, caption = caption, …) 5. pander::pandoc.table.return(…) 7. base::sapply(1:ncol(t), function(x) is.numeric(t[, x])) 8. base::lapply(X = X, FUN = FUN, …) 9. pander (local) FUN(X[[i]], …) 11. base::[.data.frame(t, , x) Warning messages: 1: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 2: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 3: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 4: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 5: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 6: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 7: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 8: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter 9: In doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a graphical parameter

pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
0.99 0.99 0.984 0.996 
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
1.07 1.07 1.07 1.07 
pander::pander(t(rrAnalysisTrain$c.index$cstatCI),caption="C. Index")
C. Index
mean.C Index median lower upper
0.777 0.777 0.758 0.796 
#pander::pander(rrAnalysisTrain$c.index,caption="C. Index")
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.824 0.804 0.845 
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.562 0.521 0.602 
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.899 0.882 0.914 
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.449 0.331 
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")

Quitting from lines 156-169 [unnamed-chunk-8] (flchain.Rmd) Error in t.default(): ! argument is not a matrix Backtrace: 1. pander::pander(t(rrAnalysisTrain\(RR_atP), caption = "Risk Ratio") 3. base::t.default(rrAnalysisTrain\)RR_atP)

pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 694.989012 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 1305 189 436.0 139.89 505.22
class=1 214 76 65.9 1.54 1.73
class=2 481 340 103.1 544.04 664.87

Performance on the test data set

index <- predict(ml,dataFLTest)
pp <- predictionStats_binary(cbind(dataFLTest$status,index),plotname="FLC")

par(op)


prob <- ppoisGzero(index,h0)
rdata <- cbind(dataFLTest$status,prob)
rrAnalysis <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataFLTest$time,
                     title="Test: FLC",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

External Data Report

pander::pander(t(rrAnalysis$OERatio),caption="O/E Ratio")

Quitting from lines 196-208 [unnamed-chunk-10] (flchain.Rmd) Error in [.data.frame: ! undefined columns selected Backtrace: 1. pander::pander(t(rrAnalysis\(OERatio), caption = "O/E Ratio") 2. pander:::pander.htest(t(rrAnalysis\)OERatio), caption = “O/E Ratio”) 3. pander::pandoc.table(res, caption = caption, …) 5. pander::pandoc.table.return(…) 7. base::sapply(1:ncol(t), function(x) is.numeric(t[, x])) 8. base::lapply(X = X, FUN = FUN, …) 9. pander (local) FUN(X[[i]], …) 11. base::[.data.frame(t, , x) There were 18 warnings (use warnings() to see them)

pander::pander(t(rrAnalysis$OE95ci),caption="O/E Ratio")
O/E Ratio
mean 50% 2.5% 97.5%
0.946 0.946 0.942 0.951 
pander::pander(t(rrAnalysis$OAcum95ci),caption="O/Acum Ratio")
O/Acum Ratio
mean 50% 2.5% 97.5%
1.05 1.05 1.05 1.05 
pander::pander(t(rrAnalysis$c.index$cstatCI),caption="C. Index")
C. Index
mean.C Index median lower upper
0.797 0.797 0.785 0.809 
#pander::pander(rrAnalysis$c.index,caption="C. Index")
pander::pander(t(rrAnalysis$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
est lower upper
0.838 0.825 0.851 
pander::pander((rrAnalysis$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
est lower upper
0.587 0.561 0.614 
pander::pander((rrAnalysis$ROCAnalysis$specificity),caption="Specificity")
Specificity
est lower upper
0.905 0.894 0.915 
pander::pander(t(rrAnalysis$thr_atP),caption="Probability Thresholds")
Probability Thresholds
90% 80%
0.449 0.331 
pander::pander(t(rrAnalysis$RR_atP),caption="Risk Ratio")

Quitting from lines 196-208 [unnamed-chunk-10] (flchain.Rmd) Error in t.default(): ! argument is not a matrix Backtrace: 1. pander::pander(t(rrAnalysis\(RR_atP), caption = "Risk Ratio") 3. base::t.default(rrAnalysis\)RR_atP)

pander::pander(rrAnalysis$surdif,caption="Logrank test")
Logrank test Chisq = 1859.709512 on 2 degrees of freedom, p = 0.000000
  N Observed Expected (O-E)^2/E (O-E)^2/V
class=0 2953 376 992 382.6 1440
class=1 472 184 138 15.3 17
class=2 1099 797 227 1433.1 1746